function logp = gmm_logpdf(X, gmm)
    N = size(X, 1);
    d = size(X, 2);
    num_components = size(gmm.mu, 1);
    
    logp = zeros(N, 1);
    
    logdets = zeros(1, num_components);
    sigmainvs = cell(num_components, 1);
    for k = 1 : num_components
        logdets(k) = logdet(gmm.sigma{k}, 'chol');
        sigmainvs{k} = inv(gmm.sigma{k});
    end
    
    component_logp = bsxfun(@plus, zeros(N, num_components), gmm.mixing_pi);
    component_logp = component_logp - d/2 * log(2 * pi);
    component_logp = bsxfun(@plus, component_logp, -1/2 * logdets);
    
    for k = 1 : num_components
        V = bsxfun(@minus, X, gmm.mu{k}');
        component_logp(:,k) = component_logp(:,k) - 0.5 * dot(V, V * sigmainvs{k}, 2);
    end
    
    logp = logsumexp(component_logp, 2);
    
    % Previous version, not vectorized
    %
    %for i = 1 : N
    %    x = X(i,:)';
    %    component_logp = log(gmm.mixing_pi);
    %    for k = 1 : num_components
    %        component_logp(k) = component_logp(k) - d/2 * log(2 * pi) -1/2 * logdets(k);
    %        x = X(i,:)';
    %        component_logp(k) = component_logp(k) - 0.5 * (x - gmm.mu{k})' * sigmainvs{k} * (x - gmm.mu{k});
    %    end
    %    logp(i) = logsumexp(component_logp);
    %end
end